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1.  INTRODUCTION 


The  parabolic  equation  method  was  first  introduced  by  Leontovich  and  Fock 
(1944)  in  their  studies  of  tropospheric  radio  wave  propagation  over  a  long 
distance.  They  were  concerned  with  calculating  the  diffraction  of  radio  waves 
by  the  spherical  shape  of  the  earth.  The  predominate  direction  of  the  wave 
propagation,  which  is  needed  to  make  the  small -angle  parabolic  approximation, 
was  the  line  of  sight  between  the  antenna  and  the  horizon.  This  method  was 
later  applied  to  many  other  radio  wave  diffraction  problems  such  as  high 
frequency  scattering  by  obstacles  of  various  shapes. 

Applications  of  the  parabolic  equation  method  were  quickly  extended  to 
wave  propagation  problems  in  other  fields  of  physical  sciences,  such  as 
nonlinear  optics  (Svelto,  1974),  plasma  physics  (Karpman,  1975),  geophysics 
(Claerbout,  1976),  and  underwater  acoustics  (Tappert,  1977).  Efficient 
numerical  algorithms  have  also  been  developed  to  solve  parabolic  equations, 
i.e.,  linear  and  nonlinear  Schrodinger  equations,  for  large  scale  problems. 

Adoption  of  the  parabolic  equation  method  in  computations  of  linear 
surface  water  wave  propagation  started  in  the  late-70's.  Based  on  the 
parabolic  approximation,  Liu  and  Mei  (1976)  presented  an  analytical  solution 
for  the  wave  field  in  the  neighborhood  of  a  breakwater  located  on  an  uniform 
sloping  beach.  Rigorous  derivations  of  the  parabolic  wave  equations 
describing  the  general  combined  refraction  and  diffraction  were  given  by 
Radder  (1979)  and  Lozano  and  Liu  (1980).  Lozano  and  Liu  employed  the  multiple 
scales  perturbation  method  and  established  the  relationship  between  ray  theory 
and  the  parabolic  method  equation.  On  the  other  hand,  using  an  operator 
splitting  technique  (Corones,  1975),  Radder  demonstrated  that  the  parabolic 
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equations  could  be  derived  from  the  mlld-slope  equation  (Berkhoff,  1972). 


The  development  of  parabolic  equation  methods  In  water  waves  has  grown 
rapidly  In  last  five  years.  Weak  reflection  (Liu  and  Tsay,  1983), 
nonlinearity  (Kirby  and  Dalrymple,  1983;  Liu  and  Tsay,  1984),  energy 
dissipation  (Liu,  1986;  Kirby  and  Dalrymple,  1986),  and  wave-current 
interactions  (Boolj,  1981;  Liu,  1983;  Dalrymple,  Kirby  and  Hwang,  1984;  Kirby, 
1984)  have  all  been  Included  in  the  formulation.  The  parabolic  equation 
methods  have  also  been  applied  to  several  field  problems  (Dingemans,  1983;  Liu 
o>  1985).  Many  publications  concerned  with  the  parabolic  approximation 

n  elated  techniques  nave  been  summarized  In  a  recent  annotated  bibliography 
(Liu,  et  al.,  1986). 

In  this  report,  a  brief  review  of  the  mlld-slope  equation  and  parabolic 
equations  suitable  for  water  waves  Is  provided.  Special  emphasis  is  given  to 
the  physical  meanings  and  justification  of  the  parabolic  approximation. 


2.  MILD-SLOPE  EQUATION 


In  this  section,  we  present  a  simple  derivation  of  the  mild-slope 
equation  describing  the  propagation  of  small  amplitude  waves  over  a  slowly 
varying  topography.  The  derivation  is  very  similar  to  that  given  by  Smith  and 
Sprinks  (1975).  The  reader  is  also  encouraged  to  consult  the  work  of  Berkhoff 
(1972)  for  a  different  derivation. 

2.1  Derivation 


We  consider  wave  propagation  over  a  gradually  varying  topography, 

2  =  -h  (x,y).  The  Cartesian  coordinates  (x,y,z)  are  fixed  on  the  undisturbed 
water  surface  and  the  free  surface  displacement  of  small  amplitude  waves  is 
described  by  z  =  r,  (x,y,t).  Assuming  that  the  fluid  is  inviscid  and 
incompressible  and  the  flow  is  irrotational ,  we  introduce  a  velocity  potential 
$>  (x,y,z,t),  so  that  the  velocity  vector,  q,  is  defined  as  the  gradient  of  the 

■f 

potenti al ;  i  .e. ,  q  =  V4>. 

Consider  small  amplitude  monochromatic  waves  with  a  radian  frequency  w. 
The  free  surface  boundary  conditions  can  be  linearized  because  of  the 
smallness  of  the  free  surface  displacement  and  the  associated  wave  motions. 
Therefore,  the  time  dependency  in  the  surface  elevation  and  the  potential  can 
be  separated  as  follows: 


$(x,y,x,t) 


( x  ,y ,  z )  e 


-iojt 


(2.1) 


-i  wt 
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dx,y,t)  =  n(x,y)e 


(2.2) 


where  only  the  real  parts  have  physical  meaning.  The  linearized  boundary 
value  problem  for  $  can  be  written  as 


+  V24>  =  0 


-h  <  z  <  0 


(2.3) 


9d>  0) 

_  -  —  =  0 
az  g 


z  =  o 


(2.4) 


—  -  -  vh  •  V6 
az 


z  =  -h 


(2.5) 


in  which  v  =  (  JL  ,  JL  )  Is  the  two-dimensional  gradient  operator  .  Once  the 
3x  3y 

velocity  potential  is  obtained,  the  free  surface  displacement  can  be  found 
from  the  following  kinematic  free  surface  boundary  condition: 


n  =  «j>  (x,y,0) 

g 


(2.6) 


For  the  special  case  of  constant  water  depth,  the  right-hand  side  of 
(2.5)  vanishes.  The  boundary  value  problem  becomes  a  homogeneous  one.  Since 
the  water  depth  is  uniform,  there  is  no  wave  refraction.  To  satisfy  the 
boundary  conditions  (2.4),  (2.5)  and  (2.6)  we  may  rewrite  the  potential  as 


♦  -■linf 


(2.7) 


where 


cosh  k ( z+h ) 
cosh  kh 


(2.8) 


and  k  is  the  wave  number  which  is  the  solution  of  the  dispersion  relation 


to2  =  gk  tanh  kh 


(2.9) 


From  the  Laplace  equation  (2.3),  the  free  surface  displacement  must  satisfy 


the  Helmholtz  equation 


v2n  +  k2  n  =  0  (2.10) 

which  is  a  two-dimensional  reduced  wave  equation  describing  wave  diffraction 
only,  since  k  is  a  constant. 

For  slowly  varying  water  depth,  (2.7),  (2.8)  and  (2.9)  are  still  valid, 
with  k  and  h  referring  to  their  local  values.  The  governing  equation  for  n 
should,  however,  be  modified  to  include  both  refraction  and  diffraction.  A 
simple  method  of  ojtaining  the  modified  equation  is  presented  herein. 

The  function  f  satisfies  the  following  set  of  equations  for  slowly 
varying  h: 

lif  -  k2  f  =  0  ,  -h  <  z  <  0  (2.11 

3Z2 

3f  _  jh!  f  =  o  ,  z  =  0  (2.12 

32  g 

If  =  0  ,  z  =  -h  (2.13 

3Z 

provided  that  the  dispersion  relation,  (2.9),  is  also  satisfied.  Considering 
(2.3)  as  an  ordinary  differential  equation  in  z,  and  applying  Green's  second 
identity  for  <j>  and  f  in  the  region  -h  <  z  <  0,  we  have 

/°  (  >!±  f  .  iff  +)dz  •  (ii  f  -  if  »>i  -  (i*  f  -  if  *)| 

-h  3Z2  3Z 2  3Z  3z  'z=o  3Z  3Z  z  =  -h 

Upon  substituting  (2.3)  -  (2.5)  and  (2.11)  -  (2.13)  into  the  above  equation, 
we  obtain 


f°  (k- 


f  +  f  v2 )  dz  =  (f  vh  •  v<<>)  I 


(2.14) 


We  now  substitute  (2.7)  and  (2.8)  into  (2.14)  and  note  that 


7$  =  -  ( f  vn  +  n  —  vh) 

co  3h 


v2 <t>  =  -  —  [f  v2n  +  2n  —  7n  •  vh  +  n  — L  ( 7h) 2  +  n  —  72h  1 

co  3h  3h 2  3h 


32f 


3f 


Equation  (2.14)  may  be  rewritten  as 


f°  (f2  V2n  +  2  f  —  7n  •  7h  +  ti  f  ( vh ) 2 
-h  3h  ah2 

3  f  9^ 

+  n  f  —  v2h  +  k2f 2n }  dz  =  [f2  7h  •  Vn  +  n  f  —  ( vh ) 2 ] 

3h  3h  1 z=-h 


Upon  using  Leibniz'  rule,  the  first  two  terms  on  the  left  of  the  above 
equation  may  be  combined  with  the  first  term  on  the  right,  yielding 


7  •  (/  f2  dz)  7ti  +  k2  (/"  f2  dz)  ri  =  f  II  I  n  ( vh): 
-h  -h  3h  1  z=-h 


f° 

-h 


3f 


o  a2f  „  af 
-  n  /  [— y  (vh)2  +  —  v2h]  f  dz 
-h  3hz  ah 


(2.15) 


Because  the  water  depth  varies  gradually  within  the  distance  of  a  typical 
wavelength,  i.e.  |vh|/kh  <<  1,  the  right  side  of  (2.15)  is  proportional  to 
Q(|vh|/kh)2,  and  can  be  neglected.  Integrating  f2  from  z  =  -h  to  z  =  0,  we 
obtai n 

v  •  (CCg  vn)  +  k2  CCg  n  =  o  (2.16) 

where 


C=^  ,  C  =  1  m  (l  +  2k  h  ) 

k  g  2  k  sinh  2kh 

are  the  phase  velocity  and  the  group  velocity,  respectively. 


(2.17) 

Equation  (2.16) 


-8- 


describes  wave  fields  in  which  both  refraction  and  diffraction  are 


considered.  Because  of  the  assumption  of  gradually  varying  topography  adopted 
in  derivation,  (2.16)  is  now  called  the  mild-slope  equation. 

In  the  limiting  case  of  arbitrary  constant  depth  including  the  deep  water 
case  where  kh  >>  1,  (2.16)  reduces  to  the  Helmholtz  equation  (2.10).  On  the 
other  hand,  in  the  shallow  water  limit,  kh  <<  1,  (2.16)  becomes 

v.(hvn)+  —  n  =  0  (2.18) 

9 

since  C  =  Cg  =  /gh  and  w2  =  gk2h.  Equation  (2.18)  is  the  linear  shallow  water 
wave  equation  and  can  be  derived  directly  from  the  continuity  and  momentum 
equations  without  using  the  mild-slope  assumption.  In  other  words,  (2.18)  is 
valid  even  for  0(|vh|/kh)  =  0(1).  Therefore,  the  mild-slope  equation,  (2.16), 
provides  an  interpolation  for  the  whole  range  of  wavelength. 

Using  the  following  transformation: 

n  =  s//CCg  (2.19) 

in  (2.16),  we  obtain 

vH  +  kc2  5  =  0  (2.20) 

where  the  effective  wave  number  kc  is  defined  by 


^c 


2 


=  k2 


v2(/CCg) 

/CCg 


(2.21) 


The  second  term  on  the  right-hand  side  of  (2.21)  is  of  the  order  of 
0( j vh | /kh ) 2 .  Therefore,  to  be  consistent  with  the  order  of  magnitude  of  the 
accuracy  of  the  mild-slope  equation,  the  effective  wave  number  can  be 
approximated  by  the  local  wave  number  k.  Thus 


V2C  +  k25  =  0 


(2.22) 


which  is  the  Helmholtz  equation  with  a  variable  coefficient  k  being  determined 
by  the  dispersion  relation. 

2.2  Relation  Between  Solutions  of  the  Mi  Id-Slope  Equation  and  Ray  Theory 

In  this  section,  we  examine  the  relationship  between  solutions  of  the 
mi  Id- si  ope  equation  and  those  of  ray  theory.  To  facilitate  the  analysis,  we 
rewrite  the  free  surface  displacement  as 

n  =  A  (x,y)  elS  (2.23) 

where  both  A  and  S  are  real  functions.  Substituting  (2.23)  into  (2.16), 
we  obtain 

i  [vS  •  (CCg  vA)  +  V  .  (CCg  AvS)]  +  7  .  (CCg  vA) 

-  vS  .  (CCg  A7S)  +  k2  CCg  A  =  0 
Multiplying  the  above  equation  by  A  yields 


VS 


k2)  A 


A  [CCg  V2A  +  7 ( CCg )  .  VA  -  CCg  (| 

+  i  v  •  (CCg  A2tS)  =  0 

The  real  and  imaginary  parts  of  (2.24)  should  be  zero.  Thus 


vS  2  -  k2 


v2A  v(CC„ )  •  vA 


CCg  A 


and 

7  •  (CCg  A2  7S)  =  0 
Let  us  define 


(2.24) 


(2.25) 


(2.26) 
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and  substitute  this  into  (2.26)  to  get 


‘  (‘.  A2)  *  ‘fsj  V  7  (1^)1  A2  ' 


(2.28) 


Equation  (2.25)  reduces  to  the  well-known  eikonal  equation  of  ray  theory,  if 
the  right-hand  side  terms  are  ignored.  At  the  same  time,  (2.28)  becomes  the 
transport  equation  which  requires  the  conservation  of  wave  action  A2/u>  along  a 
wave  ray. 

Therefore,  by  solving  the  mild-slope  equation,  the  effects  of 
diffraction,  i.e.,  the  gradient  and  the  curvature  of  the  amplitude,  are  taken 
into  account.  Wave  action  is  no  longer  conserved  along  a  "ray".  Instead, 
(2.25)  and  (2.28)  show  that  there  is  a  gain  or  loss  in  wave  action  along  a 
"ray"  because  of  diffraction. 


2.3  Energy  Dissipation 


In  the  previous  sections  the  mild-slope  equation,  (2.16),  and  the  ray 
theory  were  derived  based  on  the  assumption  that  no  energy  dissipation  occurs 
during  the  wave  propagation  process.  However,  in  most  coastal  engineering 
problems  the  energy  dissipation  effects,  such  as  bottom  friction  and  wave 
breaking,  may  become  important.  The  mild-slope  equation  can  be  modified  in  a 
simple  manner  to  accommodate  these  phenomena. 

Let  us  consider  Wr  (x,y)  as  an  energy  dissipation  function  describing  the 
rate  of  change  of  wave  energy.  The  energy  equation,  (2.28),  may  be  modified 
to  be 
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7  .  (C  A2)  +  [——  C  •  V  (-L — L)  ]  A2  -  W  A2  =  0 
g  | VS|  g  k  I  VS|  r 


(2 


Since  energy  dissipation  also  affects  the  phase  of  the  wave  train,  the 
equivalent  eikonal  equation,  (2.25),  should  also  be  modified  to  be 


VS 


2  _  k2  =  v2A  +  V(CCg)  .  VA  + 


CCg  A 


CCr 


(2 


where  Wf  is  a  prescribed  energy  dissipation  function  influencing  the  phase 
function  S.  Combining  (2.29)  and  (2.30)  and  introducing  a  complex  energy 
dissipation  function 

W  =  Wr  +  i  Wf  (2 

we  find  the  corresponding  mi ld-slope  equation 

v  •  (CC  Vn)  +  k2  CC  (1  -  1W_)  n  =  0  (2 

9  9  kCg 

The  energy  dissipation  function  is  determined  according  to  different 
dissipative  processes  and  could  be  a  function  of  the  free  surface 
displacement,  n.  which  makes  (2.32)  nonlinear  (Dalrymple,  Kirby  and  Hwang, 
1984;  Kirby  and  Dalrymple,  1986). 

Using  the  same  transformation  as  (2.19)  for  n,  we  obtain 


v2c  +  kc2  c  =  0 


(2 


with 


k 


2 


C 


V2(/CCg)  ikW 

^5  cg 


(2 


If  energy  dissipation  is  ignore i,  (2.33)  reduces  to  (2.21). 


3.  PARABOLIC  APPROXIMATION  AND  PARABOLIC  WAVE  EQUATIONS 


The  mild-slope  equations,  (2.16)  and  (2.22),  have  been  used  to  study 
various  kinds  of  combined  refraction  and  diffraction  phenomena  as  boundary 
value  problems  (e.g.,  Smith  and  Sprinks,  1975;  Skovgaard  and  Jonsson,  1980; 
Houston,  1981;  Tsay  and  Liu,  1983).  It  is  not,  however,  always  convenient  to 
treat  every  problem  as  a  boundary  value  problem,  especially  if  the  boundary 
conditions  are  not  well  defined  along  the  entire  boundary.  For  instance,  if  a 
shoreline  is  part  of  the  boundary,  the  boundary  conditions  for  breaking  waves 
along  the  shoreline  are  not  certain.  Furthermore,  if  the  mi ld-slope  equation 
is  solved  directly,  the  size  of  the  finite  elements  (or  grid  cells)  used  in 
discretizing  the  computational  domain  must  be  a  small  fraction  of  the  local 
wavelength.  A  huge  number  of  elements  (or  nodes)  may  be  necessary  when  the 
computational  domain  is  large.  An  alternative  computing  method  is  needed  to 
address  regional  (large  scale)  problems  concerning  wave  propagation  in  coastal 
waters.  The  parabolic  equation  method,  discussed  in  the  next  section, presents 
a  possible  alternative. 

3.1  Forward  Scattering 

The  basic  concept  of  the  parabolic  equation  method  is  to  convert  the 
mi ld-slope  equation  to  a  set  of  approximate  equations  which  describe  a  wave 
system  propagating  in  a  prescribed  direction  while  still  considering  the 
energy  flux  in  the  transverse  direction.  Thus,  both  refraction  and 
diffraction  are  included  in  this  approximate  formulation.  There  are  different 
ways  to  obtain  the  parabolic  equations,  such  as  the  matrix  splitting  method 


(Corones,  1975;  Radder,  1979)  and  the  multi -scale  perturbation  method  (Lozano 
and  Liu,  1980).  In  this  section,  we  present  a  simple  approach,  primarily  due 
to  Claerbout  (1976)  and  Tappert  (1977),  to  illustrate  the  method. 

We  assume  that  the  primary  wave  propagation  direction  is  in  the 
x-direction.  The  modified  free  surface  displacement  £,  (2.19),  can  be  written 
as 

5  =  y  (x,y)  eikox  (3.1) 

where  k0  is  a  reference  wave  number  and  ¥  is  a  slowly  varying  function  of  both 
x  and  y  provided  that  the  difference  between  k0  and  k  is  small.  Substituting 
(3.1)  into  (2.20),  we  obtain 


«\2m  2  jy  'syu  2  2 

U  +  11  +  2ik0  —  +  (kc  -  k0)  *  =  0 

ax2  ay2  ax 


(3.2) 


For  the  case  of  constant  water  depth  without  energy  dissipation,  k0  =  kc=  k, 
and  (3.2)  becomes 


2ik  ”  -  0 
ax 


(3.3) 


The  usual  parabolic  approximation  assumes  that  the  length  scale  of  the 
amplitude  variation  in  the  x-direction  (direction  of  wave  propagation),  Lx, 
is  much  longer  than  the  length  scale  of  the  amplitude  variation  in  the 
y-direction  (transverse  direction),  Ly.  Because  ¥  is  a  slowly  varying 
^unction  in  x  and  y,  the  length  scales  Lx  and  Ly  are  much  longer  than  the 
typical  wavelength.  In  other  words, 

0  =  e  «  1  , 

Lx 


-14- 


Hr 


0  (k  L  )  =  0  Ul)2  =  0  (e  ‘) 

X  Ly 


0  (kL  )  =  0  (—'  =  0  (e"1) 


(3.4) 


Because  of  these  relationships  among  the  length  scales,  32f/3x2  is  much 


smaller  than  32,F/3y2  and  can  be  neglected.  Therefore,  (3.2)  can  be 


approximated  as 


+  2i k 0  —  +  (kc  -  kn)  y  =  0 


(3.5) 


The  order  of  magnitude  of  the  term  dropped  from  (3.5)  is  of  0(  e2).  If  the 


third  term  in  (3.5)  is  in  the  same  order  of  magnitude  as  the  other  two  terms, 


it  is  required  from  (3.4)  that 


(-i)2  -  1  =  0( e2 ) 
ko 


(3.6) 


If  the  water  depth  is  a  constant  and  energy  dissipation  is  ignored,  (3.3)  can 


be  simplified  farther  to  be 


if*  +  2ik  2!  =  0 


(3.7) 


the  above  equation  is  in  the  same  form  as  the  heat  equation  if  x  is 


interpreted  as  time  and  the  heat  conductivity  is  taken  to  be  imaginary. 


Equations  (3.5)  and  (3.7)  have  the  form  of  the  linear  Schrodinger  equation. 


To  examine  further  the  approximations  involved  in  the  parabolic  wave 


equation,  we  can  convert  (3.5)  or  (3.7)  back  from  the  v  variable  to  the  £ 


variable  and  compare  them  with  the  Helmholtz  equation,  (2.22);  i.e., 


v2C  +  k2e  =  0 


(2.22) 


For  simplicity,  we  focus  our  discussion  on  the  constant  water  depth  without 
energy  dissipation  case.  Substitution  of  (3.1)  into  (3.7)  yields 

!i  +  2ik  ii  +  2k 2 c  =  0  (3.8) 

3y2  3x 

which  is  significantly  different  from  the  Helmholtz  equation.  The  possible 
solutions  for  the  above  equation  can  be  expressed  as  exp[i(kx  x  +  ky  y ) ] , 
which  represents  a  plane  wave  propagating  in  the  k  =  (kx,  ky)  direction. 
Substituting  the  exponential  function  solution  into  (3.8),  we  obtain 
an  algebraic  equation  for  kx  and  ky: 

ky2  +  2  k  kx  -  2  k2  =  0  (3.9) 

Note  that  kx  is  not  equal  to  k  unless  ky  is  zero.  The  difference  between  k 
and  kx  is  partly  accounted  for  by  the  phase  modification  in  Y,  (3.1). 

Equation  (3.9)  is  plotted  in  Figure  3.1,  which  shows  a  parabolic  curve.  Since 
only  positive  kx  values  are  allowable,  the  parabolic  equation,  (3.7), 
describes  only  forward  scattering. 

Substituting  the  exponential  solution,  exp  [i(kx  x  +  ky  y)  ]  into  the 
Helmholtz  equation,  we  obtain  the  relation  between  ky,  kx  and  k  as 

kx2  +  ky2  =  k2  (3.10) 

which  represents  a  circle  on  the  kx ,  ky  plane  (Figure  3.1).  It  is  clear  from 

the  figure  that  (3.9)  is  an  approximation  of  (3.10)  for  small  angle,  8  = 

-1 

tan  (ky/kx).  Furthermore,  (3.10)  indicates  that  for  a  given  ky  the  sign  of 
kx  could  be  either  positive  or  negative.  In  other  words,  from  the  Helmholtz 
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equation  waves  may  propagate  in  either  the  positive  or  negative  x-di recti  on 
with  the  same  wave  number  kx;  whereas,  the  parabolic  approximation,  (3.7), 
examines  only  forward  scattering,  with  an  additional  small  angle  requirement 
The  small  angle  requirement  may  be  relaxed  if  the  forward  propagating 
wave  equation  has  a  semi -circle  for  its  relation  between  kx  and  ky.  The 
equation  for  a  perfect  semi -circle  is  given  by 

kx  =  /kz  -ky2  (3.11 

The  corresponding  wave  equation  is 

i  i£  +  k  0  5  =  0  (3.12 

ax 

where  Q  is  a  pseudo-differential  operator  given  by 

1  a2  1/2 

Q  =  (l  +  -L  *  )  '  (3.13 

kz  ay 2 

(One  can  easily  obtain  (3.11)  by  substituting  the  exponential  solution, 
exp  [i(kx  x  +  ky  y)]  into  (3.12)).  As  indicated  in  (3.4),  for  the  parabolic 

approximation  where  0(kLy)  -  OU"1),  (3.13)  can  be  written  as 

Q  =  (1  +  e2)1/2  (3.14 

For  small  e,  0  may  be  approximated  using  several  different  methods.  For 
instance,  using  a  binomial  expansion,  for  small  e,  0  may  be  written  as 

0  =  (1  +  _LlL  +  +  ...)  (3.15 

2k2  ay2  8k4  ay4 

The  corresponding  wave  equation  becomes 

i  15.  +  k  (1  +  J_  iL  +  _L  i!_  +  ...)  c  =  0  (3.16 

ax  2k2  ay2  8k4  ay4 
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By  substituting  the  exponential  solution,  exp  [i ( kx  x  +  ky  y)  ]  into  (3.16), 
the  associated  equation  for  kx  and  ky  may  be  written  as 


k  (1 


k  2 

Ky 

2k2 


kv 

-  ...) 

8k4 


(3.17) 


which  is  the  binomial  expansion  of  (3.11)  for  a  small  angle  tan  (ky/kx). 
This  expansion  converges  for  all  0  <  ky  <  k.  Therefore,  the  parabolic 
approximation  shown  in  (3.8)  and  (3.9)  is  the  lowest  order  approximation  of 
(3.16)  and  (3.17),  being  truncated  at  the  second  term  of  the  series 
expansion.  The  effects  of  large  angle  diffraction  may  be  studied  if  the 
higher  order  terms  are  kept  in  the  analysis.  In  table  3.1,  a  comparison 
between  the  exact  solution,  (3.11),  and  the  approximated  solutions  with  two 
terms  and  three  terms  is  given. 

As  an  alternative  to  the  binomial  expansion,  we  can  choose  the  rational 
expansion  of  Q,  for  small  e, 


0  = 


(1  +  e2 ) 


V2 


=  (1  + 


ay' 


r)/(l  + 


(3.18) 


Substituting  (3.18)  into  (3.12)  yeilds  the  following  wave  equation 


(  ii  +  *  kC .  Jiis  =  o  (3.i9) 

ax  4k2  axay2  4k  ay2 


The  correspond' ng  equation  for  kx  and  ky  can  be  written  as 
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The  accuracy  of  (3.20)  compared  with  the  exact  solution  is  also  shown  in  table 
3.1.  It  is  clear  that  the  approximation  given  in  (3.20)  is  the  most  accurate 
one  for  larger  angles. 


Substituting  (3.1)  into  (3.16)  and  (3.19)  with  k=k0,  the  higher-order 
parabolic  approximations  can  be  written  in  terms  of  y  as 


2ik 

SY  + 

a2'!' 

-  1  -  0 

(3.21) 

3X 

3y2 

% 

Oi 

2ik 

3Y 

32f 

«■  1  s3t  -  0 

(3.22) 

3X 

ay2 

2k  3x ay2 

respectively.  In  fact,  (3.21)  can  be  obtained  from  (3.22)  by  using  the  first 
order  approximation,  (3.7),  of  3 ^/ax  in  the  third  term.  If  the  water  depth  is 
a  function  of  spatial  coordinates  or  the  energy  dissipation  is  not  negligible, 
kc  is  no  longer  a  constant.  The  higher-order  parabolic  wave  equation  can  be 
written  as  follows: 


2ikn  1!  ♦ 


32¥ 


+  (k l  -  kn2)  *  + 


a3f 


3x  ay* 


2k  0  3x  3y  2 


(3.23) 


Note  that  if  the  last  term  in  the  above  equation  is  ignored,  the  lower-order 
parabolic  approximation  equation,  (3.5),  is  recovered.  Using  a  different 
approach  as  shown  in  the  next  section,  the  parabolic  wave  equations  can  be 
expressed  in  slightly  different  forms: 
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f  -  . 
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i  a3* 
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0 


(3.24) 


We  remark  here  that  since  kc  -  k0  =  0{e2)  from  (3.6),  (3.24)  can  be  reduced  to 
(3.23)  by  replacing  kc  with  k0  without  losing  accuracy. 

The  lower-order  (small  angle)  parabolic  wave  equations,  (3.5)  and  (3.7), 
are  in  the  form  of  heat  equation.  Therefore,  the  wave  propagation  problem  can 
now  be  viewed  as  an  initial  boundary  value  problem  treating  the  direction  of 
wave  propagation  as  the  time.  The  initial  conditions  in  terms  of  wave 
amplitudes  and  phase  functions  must  be  prescribed  along  a  straight  line  x  = 
constant.  The  Crank-Nicol son  finite  difference  scheme  (e.g.,  Smith,  1978)  can 
be  used  to  discretize  the  parabolic  wave  equations.  This  scheme  solves  for 
the  unknowns,  y,  along  a  line  normal  to  the  direction  of  wave  propagation  one 
step  at  a  time.  A  detailed  documentation  on  this  numerical  scheme  can  be 
found  in  Tsay  and  Liu  (1986).  Two  lateral  boundary  conditions  are  required  at 
the  end  points  of  the  line.  Two  typical  types  of  boundary  usually  occur  in 
the  combined  refraction  and  diffraction  problems:  (1)  a  rigid  boundary  and 
(2)  an  open  boundary.  In  the  case  of  a  rigid  boundary,  y  =  yfo(x),  the  no-flux 
boundary  condition  requires 

n  •  Vn  =  0  on  y  =  y^ ( x )  (3.25) 

where  n  is  the  unit  normal  along  the  solid  boundary.  Using  (2.19)  and  (3.1) 
in  (3.25)  and  employing  the  n  =  v[y-yb (x ) ] / | v[y-yb (x ) ] I ,  we  obtain 


a? 

ay 


at 

[ —  +  ikf 
ax 


f  acc  y  acc 

- 2]  - _ 

2  ax  2  ay 


9  =  0 


y  =  y.  (x) 

b 


(3.26) 
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Further  simplification  may  be  made  depending  on  the  slope  of  the  solid 
boundary.  In  the  case  of  an  open  boundary,  if  the  boundary  is  located 
sufficiently  far  away  from  the  region  of  interest  and  the  topography  is 
assumed  to  be  uniform  in  the  y-direction  near  the  open  boundary,  the 
approximate  boundary  condition  can  be  written  as 


3y 

—  =  0 

ay 


y  =  y0 


(3.27) 


The  question  concerning  accurate  and  efficient  numerical  method  for  dealing 
with  the  open  boundary  condition  remains  open.  The  reader  is  encourged  to 
consult  the  work  of  Cohn  and  Jennings  (1983),  which  discussed  different 
methods  for  treating  the  silent  open  boundary. 

Because  f  is  a  complex  quantity,  it  may  be  written  as 


y  =  Me 


-1  Im  v 

g  =  tan  ( - ) 

Re  y 


The  resulting  wave  number  components  in  the  x-  and  y-direction  can  be 
expressed  as  (k0  +  ae/ax,  ae/ay).  The  directions  of  wave  propagation  can  be 
calculated  accordingly. 


3.2  Weak  Reflection 


In  the  previous  section  only  the  forward  propagating  wave  field  was 
considered.  If  wave  reflection  is  important,  several  different  approacheshave 
been  developed  to  derive  a  set  of  parabolic  equations  for  both  forward  and 
backward  scattering  wave  fields.  The  method  presented  here  is  similar  to  that 
given  by  Tappert  (1977). 

We  first  rewrite  the  Helmholtz  equation  in  the  following  form: 


where 


I 


rtf’**'*  *  >  *-«.  j\  •%  I 


I 


_  +  kc2  Q2c  =  0 


(3.28) 


.  1  32  , 

Q2  =  _ _  —  +  1 


(3.29) 


in  which  kc  is  given  in  (2.33).  We  remark  here  that  in  the  case  of  a  constant 
water  depth  and  zero  energy  dissipation,  kc  =  k,  and  reflection  can  still  be 
caused  by  the  appearance  of  surface  piercing  structures.  Equation  (3.29) 
reduces  to  (3.13)  in  the  limiting  case. 

We  assume  that  the  total  wave  field  can  be  split  into  a  transmitted 
(forward  scattering)  wave  field  and  a  reflected  (backward  scattering)  wave 
field.  Thus 


C  =  5+  + 


—  =  ik  Q  U  -  C  ) 

3X  C  + 


From  the  equations  above,  we  obtain 


1  U  -  i_ii) 

2  kcQ  3x 


(3.30) 


(3.31) 


1  U  +  !_!£) 

2  kcQ  3x 


(3.32) 


(3.33) 


Differentiating  (3.32)  and  (3.33)  with  respect  to  x  and  using  the  governing 
equation  (3.28)  for  32?/3x2  and  (3.31)  for  3£/3x,  we  obtain 


1  3kcQ  i  a 

-  ikcQc+  + - —  c+  = - —  c- 


1  3kcQ 


2kcQ  3  x 


2kc0  3  x 


(3.34) 


-  +  ikc05_  + 


(3.35) 


2kcQ  3  x  2kcQ  3  x 


Using  the  rational  expansion  of  the  pseudo-differential  operator,  Q,  (3.29 
in  (3.34)  and  (3.35),  we  can  derive  the  following  approximate  equations: 


H+  1  335+  3i  32c+  1  3kc  1  3kc 

-  -  ikcc+  + - - - + - 5+  = - 5_  (3.36) 

3x  4kc2  3x3 y2  4kc  3y2  2kc  3x  2kc  3x 


35.  1  335_  3i  325.  1  3kr  1  3kr 

— -  +  ikc5_  + - _  + - —  + - £  5_  = - £  5+  (3.37) 

3x  4kc2  3x3y2  4kc  3y2  2kc  3x  2kc  sx 


We  remark  that  the  left-hand  side  of  the  governing  equation  for  the  forward 
scattering  wave  field  is  exactly  the  same  as  that  of  (3.19)  for  a  constant 
water  depth  and  no  energy  dissipation  (kc  =  k).  Equations  (3.36)  and  (3.37) 
denote  a  set  of  coupled  equations  for  forward  and  backward  scattering  wave 
fields.  In  terms  of  the  amplitude  functions 


*+  -  C, 


+i  k  0x 


(3.38) 


equations  (3.36)  and  (3.37)  can  be  rewritten  as 


3kc  i  3J¥+ 

+  [2kc  (kc  -  k0)  +  i  - ]  ?  + - - 


(3.40) 


3  k0  ,  32y_  .  3kc  2iknx 

-  (_  -  _2_)  -  =  i  — £  v+  e^1K  °x 

2  2kc  3y2  3x 

In  principle,  (3.39)  and  (3.40)  can  be  solved  simultaneously  for  y+  and  f_. 
Liu  and  Tsay  (1983)  developed  an  iterative  nunerical  scheme  for  (3.39)  and 
(3.40)  with  successive  corrections  on  the  reflected  and  transmitted  wave 
fields.  Their  scheme  can  be  briefly  described  as  follows:  First  the 
reflection  is  entirely  ignored;  i.e.,  ¥  0  =  0,  and  (3.39)  becomes 


3kCn  ...  o  ,  1  9  V 


[2kc  (kc  -  k0)  +  i  _±]  ¥+°  + 


2kc  3xay2 


3  k0  324'  0 

+  (-  -  - )  - —  =  0 

2  2kc  ay2 


(3.41) 


which  is  exactly  the  same  as  (3.24),  and  can  be  solved  numerically.  Equation 
(3.41)  can  be  rewritten  as  P  [t^0]  =  0  where  P  is  the  operator  of  (3.41). 

Once  0  is  found,  we  can  substitute  it  into  the  right-hand  side  of  (3.40) 
to  get  the  modified  reflected  wave  field,  y_1.  Thus 

*  3kc  2iknx 

P  [f  M  =  i  -  v  0  e  (3.42) 

-  3x  + 

★ 

where  P  represents  the  operator  on  the  left-hand  side  of  (3.40).  We  can  use 
the  same  numerical  scheme  as  described  in  the  previous  section  to  solve  (3.42) 
for  'F_1,  which  is  in  turn  substituted  into  the  right-hand  side  of  (3.39)  to 
find  the  improved  v+x.  The  procedure  should  be  repeated  until  the  converged 


solutions  are  obtained. 


4.  WEAKLY  NONLINEAR  WAVES 


The  parabolic  approximation  can  also  be  applied  to  weakly  nonlinear 
waves.  Using  multiple  scale  perturbation  schemes,  several  researchers  (Yue 
and  Mei ,  1980;  Kirby  and  Dalrymple,  1983;  and  Liu  and  Tsay,  1984)  have  derived 
the  cubic  Schrodinger  equation  as  the  evolution  equation  describing  the  wave 
envelope  of  Stokes  second-order  waves.  Here  we  present  a  simpler  approach  to 
obtain  the  same  evolution  equation. 

First,  we  rewrite  the  Helmholtz  equation,  (2.22),  in  the  following  form: 

72£  +  m?k -  5  =  0  (4.1) 

g  tanh  kh 

The  wave  frequency  in  nonlinear  wave  theory  depends  not  only  on  wave  number, 
k,  but  also  on  the  amplitude,  a.  The  nonlinear  dispersion  relation  has  the 
general  form 

O,2  =  u)2 (k 2 ,  a2)  (4.2) 

On  the  assumption  that  the  wave  amplitude  is  fairly  small,  we  can  approximate 
(4.2)  in  a  Taylor's  series  as 

=  to 2 ( k 2 )  +  [1m2—]  a2  (4.3) 

o  3(a2)  o 

where  u>o(k2)  represents  the  linear  dispersion  relation  given  in  (2.9).  For 
Stokes  waves  [3u>2/3(a2) ]0  can  be  expressed  as  (Mei,  1983): 


3gj2  gk3  tanh  kh 
- ]  =  1 - 

3(a2)  o  8  sinh4  kh 


[8  +  cosh  4kh  -  2  tanh2  kh] 


(4.4) 
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We  also  point  out  that,  from  (2.19)  and  (3.1),  the  wave  amplitude  can  be 
expressed  in  terms  of  5  and  4*  as 

a  =  |n|  =  |5|//CCj  =  |f|//CC^  (4.5) 

Following  the  analysis  presented  in  section  3.1,  we  can  readily  obtain  the 
parabolic  equation  for  the  nonlinear  Stokes  waves  as 

—  +  2ik0  ~  +  K  1 j  2  f  =  0  (4.6) 

ay2  ax  11 

where 

K  =  - - — , -  [8  +  cosh  4kh  -  2  tanh2  khl  (4.7) 

8Cg  sinh4  kh  L  J 

Equations  (4.6)  and  (4.7)  are  identical  to  those  derived  by  Yue  and  Mei 
(1980).  The  Crank-Nicol son  scheme  can  still  be  used  to  solve  the  nonlinear 
cubic  Schrodinger  equation,  (4.6).  This  equation  has  been  used  to  examine  the 
wave  field  in  the  neighborhood  of  caustics  (e.g.,  Kirby  and  Dalrymple,  1983; 


5.  CONCLUDING  REMARKS 


In  this  report,  a  brief  review  of  the  mi  Id-slope  equation  and  the  para¬ 
bolic  wave  equations  suitable  for  water  waves  was  given.  The  physical  mean¬ 
ings  and  the  limitations  of  various  parabolic  approximations  were  discussed. 
Parabolic  wave  equations  including  the  effects  of  energy  dissipations  and  non 
linearity  were  derived. 

Recently,  Kirby  (1986)  presented  a  wave  equation  which  includes  the 
effects  of  the  small -amplitude  fast  undulations  in  seafloor.  An  additional 
term  was  added  to  the  mi ld-slope  equation.  The  parabolic  wave  equations, 
similar  to  those  presented  in  this  report,  can  be  derived. 

The  concept  of  the  parabolic  approximation  presented  in  this  report  has 
also  been  extended  to  shallow-water  wave  problems  (Liu,  Yoon  and  Kirby, 

1985).  Instead  of  the  mi ld-slope  equation,  the  Boussinesq  equation  was  used 
as  the  governing  equation.  Expressing  the  solution  in  a  Fourier  series,  Liu, 
Yoon  and  Kirby  derived  a  set  of  coupled  parabolic  wave  equations  for  Fourier 
components. 

In  a  companion  report  (Tsay  and  Liu,  1986),  various  numerical  schemes  fo 
solving  the  parabolic  equation  are  presented.  The  accuracy  and  the  lim¬ 
itations  of  each  scheme  are  discussed  and  compared. 
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